#include "mcmcRegionDynamics.h"
#include "mcmcException.h"
#include "mcEdgeOp.h"
//#include "CgEdgePrcs.h"

////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////
ModelAndEnergy::ModelAndEnergy()
{
	pmodel = NULL;
	energy_shape = 0.0;
}

ModelAndEnergy::~ModelAndEnergy()
{
	Free();
}

void ModelAndEnergy::Free(void)
{
	if (pmodel != NULL)
		delete pmodel;
	pmodel = NULL;
	energy_shape = 0.0;
}

void ModelAndEnergy::CopyFromRegion(MCMC_RegionElement *pregion)
{
	Free();

	pmodel = pregion->Model()->CopyTo();
	energy_shape = pregion->EnergyShape();
}


///////////////////////////////////////////////////
// the dynamic of splitting simple region
MCMC_Dynamic::MCMC_Dynamic(MCMC_Map &map, MCMC_FilterResp &respFilter1,
						   MCMC_ModelSupervisor &model_supervisor, int &step)
{
	pmap = &map;
	prespFilter = &respFilter1;
	num_region_layer.SetDimension(3,1);
	mx_proposal.SetDimension(prespFilter->rows(),prespFilter->cols());
	pstep = &step;
	pmodel_supervisor = &model_supervisor;
	mask_before.SetDimension(prespFilter->rows(),prespFilter->cols());
	mask_before.InitValue(0);
}

MCMC_Dynamic::~MCMC_Dynamic()
{
	pmap = NULL;
	prespFilter = NULL;
}

void MCMC_Dynamic::ChangeMapPointer(MCMC_Map &map)
{
	pmap = &map;
}



void MCMC_Dynamic::SetLayerPdf(RF_uniform &rf)
{
	double		rfp[3]={1.0,5.0,10.0},dtemp;
	int			i,num_min=99999;

	// get the minimal number of regions
	for (i=0; i<3; i++)
	{
		if (num_region_layer(i)<num_min && num_region_layer(i)!=0)
			num_min = num_region_layer(i);
	}
	// set probability
	if (num_min<99999)
	{
		for (i=0; i<3; i++)
		{
			if (num_region_layer(i)>0)
			{
				dtemp = min(rfp[i],((double)num_region_layer(i))/((double)num_min));
				rf.SetOneBin(i,1/dtemp);
			}
			else
				rf.SetOneBin(i,0.0);
		}
		rf.normalize();
	}
}

void MCMC_Dynamic::Propose(void)
{
	// set proposal probabilites to be zero
	Gw1_dw2 = 0.0;
	Gw2_dw1 = 0.0;
	rate = 0.0;
	log_pos1 = 0.0;
	log_pos2 = 0.0;
	bValidProposal = false;

	for (int i = 0; i < 4; i++){
		sprintf(Info[i], "");
	}
}

// rand generator to accept a model under certain probability
bool MCMC_Dynamic::TryAccept(void)
{
	char	pstemp1[200],pstemp2[200];
	Rand	rf_accept;

	bAccept = false;
	if (Gw2_dw1>0.0 && Gw1_dw2>0.0)
	{
		rate = Gw2_dw1/Gw1_dw2*exp(-(log_pos2-log_pos1));
		bAccept = rf_accept.Accept(min(1.0,rate));
	}

	DoubleToString(pstemp1,Gw1_dw2);
	DoubleToString(pstemp2,log_pos2);
	sprintf(Info[1], "Gw1_dw2=%s, E2=%s", pstemp1,pstemp2);
	DoubleToString(pstemp1,Gw2_dw1);
	DoubleToString(pstemp2,log_pos1);
	sprintf(Info[2], "Gw2_dw1=%s, E1=%s", pstemp1,pstemp2);
	DoubleToString(pstemp1,rate);
	sprintf(Info[3], "rate=%s, Accept=%d",pstemp1,bAccept);

	return bAccept;
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
MCMC_Dn_Switch::MCMC_Dn_Switch(MCMC_Map &map, MCMC_FilterResp &respFilter1,
							 MCMC_ModelSupervisor &model_supervisor, int &step)
:MCMC_Dynamic(map,respFilter1,model_supervisor,step)
{
	vc_order_indexid.SetRowVector(LayerRegion()->MaxGRegionIndex());
	vc_order_indexid.InitValue(-1);
}

MCMC_Dn_Switch::~MCMC_Dn_Switch()
{
	FreeModels();
}

void MCMC_Dn_Switch::FreeModels(void)
{
	pmodel_old.Free();
	pmodel_new.Free();
}

bool MCMC_Dn_Switch::PrepareToSwitchModel(void)
{
	LinkList<int>	list_indexid;

	// get all the top node regions
	LayerRegion()->GetAllTopNodeRegions(list_indexid);
	list_indexid.ToVector(vc_indexid);
	if (vc_indexid.size() <= 0)
		return false;

	// set the random variable
	MCMC_RegionElement *pregion;
	rf_indexid.Set(0, vc_indexid.size()-1, 1);
	for (int i=0; i<vc_indexid.size(); i++)
	{
		pregion = LayerRegion()->GetRegionFromIndexId(vc_indexid(i));
		assert(pregion != NULL);

		rf_indexid.SetOneBin(i, pregion->GetCryingVoice());
	}
	rf_indexid.normalize();

	for (int i=0; i<vc_indexid.size(); i++)
		vc_order_indexid(vc_indexid(i)) = i;
	return true;
}

double MCMC_Dn_Switch::ProposeRegionToSwitch(void)
{
	if (!PrepareToSwitchModel())
	{
		indexid_region = -1;
		return 0.0;
	}
	else
	{
		int i;

		i = rf_indexid.SamplingIntNonSetCdf();
		indexid_region = vc_indexid(i);
		return rf_indexid.pdf(i);
	}
}

double MCMC_Dn_Switch::Q_ProposeRegionToSwitch(void)
{
	if (!PrepareToSwitchModel() || indexid_region<0)
		return 0.0;
	else
	{
		int i;

		i = vc_order_indexid(indexid_region);
		return rf_indexid.pdf(i);
	}
}
void MCMC_Dn_Switch::Propose(void)
{
	double			q1=1.0,q2=1.0;
	bool			bChange = false;
	MCMC_RegionElement *pregion;

	// set proposal probabilites to be zero
	MCMC_Dynamic::Propose();

	FreeModels();
	q1 *= ProposeRegionToSwitch();

	pregion = LayerRegion()->GetRegionFromIndexId(indexid_region);
	if (pregion != NULL)
	{
		bValidProposal = true;

		pmodel_old.CopyFromRegion(pregion);
		// compute the likelyhood
		log_pos1 = pmodel_old.TotalEnergy()/LayerRegion()->GetT();

		LinkList<int> list_aregions;

		pregion->GetAllAtomicRegions(list_aregions);

		pmodel_new.pmodel = pmodel_supervisor->ProposeModel(list_aregions);

		pmodel_new.energy_shape += pregion->EnergyShape();
		
		// add shape prior
		log_pos2 = pmodel_new.TotalEnergy()/LayerRegion()->GetT();

		// compute the probability to switch
		DoSwitch();

		q2 = Q_ProposeRegionToSwitch();

		Gw1_dw2 = q1*pmodel_supervisor->GetPdf(pmodel_new.pmodel);
		Gw2_dw1 = q2*pmodel_supervisor->GetPdf(pmodel_old.pmodel);
	}

	sprintf (Info[0],"iteration=%d T=%.2f, switch region %d", Step(), pmap->GetT(),indexid_region);
}

void MCMC_Dn_Switch::DoSwitch(void)
{
	MCMC_RegionElement *pregion;

	pregion = LayerRegion()->GetRegionFromIndexId(indexid_region);
	if (pregion != NULL)
		pregion->SetStat(pmodel_new.pmodel);
}

void MCMC_Dn_Switch::RecoverSwitch(void)
{
	MCMC_RegionElement *pregion;

	pregion = LayerRegion()->GetRegionFromIndexId(indexid_region);
	if (pregion != NULL)
		pregion->SetStat(pmodel_old.pmodel);
}

void MCMC_Dn_Switch::DoAction(void)
{
	//cout << "Switch" << endl;
	if (bValidProposal)
	{
		if (!bAccept)
			RecoverSwitch();
	}
	for (int i=0; i<vc_indexid.size(); i++)
		vc_order_indexid(vc_indexid(i)) = -1;
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
MCMC_Dn_Diffusion::MCMC_Dn_Diffusion(MCMC_Map &map, MCMC_FilterResp &respFilter1,
							 MCMC_Diffusion &mcmcDiffusion, MCMC_ModelSupervisor &model_supervisor, int &step)
:MCMC_Dynamic(map,respFilter1,model_supervisor,step)
{
	pmcmc_diffusion = &mcmcDiffusion;
}

MCMC_Dn_Diffusion::~MCMC_Dn_Diffusion()
{
}

void MCMC_Dn_Diffusion::Propose(void)
{

	// set proposal probabilites to be zero
	MCMC_Dynamic::Propose();

	list_indexid.Empty();
	
//	list_indexid.Insert(2012);
//	list_indexid.Insert(2003);

	LayerRegion()->GetAllTopNodeRegions(list_indexid);

	layer_region_buf = *LayerRegion();

	bValidProposal = true;

	sprintf (Info[0],"iteration=%d T=%.2f, Diffusion", Step(),pmap->GetT());
}

bool MCMC_Dn_Diffusion::TryAccept(void)
{
	if (list_indexid.Num()>0)
		bAccept = true;
	else
		bAccept = false;

	return bAccept;
}

void MCMC_Dn_Diffusion::DoAction(void)
{
	if (bValidProposal && bAccept)
	{
		// diffuse the selected boundaries
		// cout << "Diffusion" << endl; // Siyuan
		try
		{
			pmcmc_diffusion->GridGo(list_indexid);

			// diffuse model parameters
			LayerRegion()->MergeSmallAtomicRegion(3);
			LayerRegion()->CreateTopologyRecords();
			LayerRegion()->CleanTopologyForGRegions();

			LayerRegion()->CreateAffinityMap(*prespFilter);
			
			pmodel_supervisor->FitAllRegionParms();
			LayerRegion()->ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);
		}
		catch (XPT_OutOfBound *pxpt)
		{
			cout << "DoAction" << endl;
			delete pxpt;

			*LayerRegion() = layer_region_buf;
		}
	}
}


/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
MCMC_Dn_RepartitionARegions::MCMC_Dn_RepartitionARegions(MCMC_Map &map, MCMC_FilterResp &resp_filter1, MCMC_ModelSupervisor &model_supervisor, mcmcMatrix<int> *partition_maps1)
{
	pmap = &map;
	pmodel_supervisor = &model_supervisor;
	presp_filter = &resp_filter1;

	for (int i=0; i<3; i++)
		p_parition_maps_edge[i] = partition_maps1[i];
/*
	MCMC_EdgeOp edge_canny;

	edge_canny = presp_filter->ras_intensity;

	Raster<double> ras_edges;
	Raster<uchar> ras_aa;

	ras_edges = edge_canny.Canny(3.0);

	ras_aa = ras_edges.ToVisibleRaster();
	SaveRaster2ImageFile(ras_aa, "a.tif");
*/	
/*
	Raster<int> ras_patition_new;
	Raster<uchar> ras_aa;

	ObtainNewPartitionMap(ras_patition_new);
	
	ras_aa = ras_patition_new.ToVisibleRaster();
	SaveRaster2ImageFile(ras_aa, "a.tif");
*/
}

MCMC_Dn_RepartitionARegions::~MCMC_Dn_RepartitionARegions()
{
}

/* Jiadi Yang
void MCMC_Dn_RepartitionARegions::ObtainNewPartitionMap(Raster<int> &ras_patition_new)
{
	mcmcMatrix<int>	mx_region_labels;
	Raster_Color<uchar> *pimage_color;
	Raster<int> ras_existing;

	pimage_color = &presp_filter->image_original;

	ras_existing.SetDimension(pimage_color->rows(), pimage_color->cols());
	ras_existing.InitValue(0);
	LayerRegion()->DrawRegions(ras_existing);

	mx_region_labels.SetDimension(pimage_color->rows(), pimage_color->cols());
	
	int num_regions;
	CEdgePrcs edge_linking_proc;	

	num_regions = edge_linking_proc.LabelRegions(3.0, mx_region_labels,true,pimage_color);

	ras_patition_new = mx_region_labels;
}
*/

void MCMC_Dn_RepartitionARegions::Propose(void)
{
	layer_region_buf = *LayerRegion();


	try
	{
		repartition.RePartitionAtomicRegions(*LayerRegion(), p_parition_maps_edge[2]);


		LayerRegion()->CreateAffinityMap(*presp_filter);
		pmodel_supervisor->FitAllRegionParms();
		LayerRegion()->ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);
	}
	catch (XPT_OutOfBound *pxpt)
	{
		cout << "Propose" << endl;
		delete pxpt;

		*LayerRegion() = layer_region_buf;
	}

/*	
	LinkList<int> list_change,list_new;
	int temp;

	temp = 0;
	list_change.Insert(temp);

	temp = 108;
	list_change.Insert(temp);

	LayerRegion()->DeleteAffinity(list_change);
	repartition.RePartitionAtomicRegions(list_new, *LayerRegion(), list_change, p_parition_maps_edge[0]);

	LayerRegion()->CreateAffinityMap(*presp_filter, &list_new);

			
	pmodel_supervisor->FitAllRegionParms();
	LayerRegion()->ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);
*/
}

